********************************************************************************

* State level data 

	foreach occ in physicians lawyers {
			
		*All ages
		
		use if year==2017 using "${clean_data}/panel_`occ'_clean.dta", clear
		format wkh %9.0g	
		
		if "`occ'" == "physicians" {
		gen primary_care=(spec_category_code==1)
		gcollapse (mean) age female married ranked_school top5_school ptotinc logptotinc totw2comp aginc logaginc top1 with25Kbizinc retired_1099ssa wkh income_acs primary_care (count) n_obs_2017_ur=age n_obs_schl_ur=ranked_school n_obs_acs_ur=income_acs, by(stfips) 
		}
		else if "`occ'" == "lawyers" {
		gcollapse (mean) age female married ptotinc logptotinc aginc wkh (count) n_obs_2017_ur=age n_obs_acs_ur=ptotinc, by(stfips) 
		}	
		ds stfips *ur, not
		foreach var in `r(varlist)' {				
				  rename `var' mean_`var'
			     }
		gen sample="`occ'; year 2017; all ages"
		order sample
		save ${temp}/state_`occ'_allage2017.dta, replace 
		
		* Ages 40 to 55 
					
		use if year==2017 & inrange(age,40,55) using "${clean_data}/panel_`occ'_clean.dta", clear
		
		if "`occ'" == "physicians" {
		merge 1:1 id year using "${output}/twfe/geo_twfe/`occ'/IndividFixedEffects_ALL_`occ'.dta", keep(master match) nogen
		}
		else if "`occ'" == "lawyers" {
		merge 1:1 id year using "${output}/twfe/geo_twfe/`occ'/IndividFixedEffects_`occ'.dta", keep(master match) nogen
		}		
		format wkh %9.0g
		gen logptotinctwfe=logptotinc if !missing(__hdfeINDIVID__)
		
		if "`occ'" == "physicians" {
		gen primary_care=(spec_category_code==1)
		gcollapse (mean) age female married ranked_school top5_school ptotinc logptotinc totw2comp aginc logaginc top1 with25Kbizinc retired_1099ssa wkh income_acs primary_care fe_cz=__hdfePLACE__ fe_id=__hdfeINDIVID__ logptotinctwfe (count) n_obs_2017_ur=age n_obs_schl_ur=ranked_school n_obs_acs_ur=income_acs n_obs_twfe_ur=__hdfeINDIVID__, by(stfips) 
		}
		else if "`occ'" == "lawyers" {
		gcollapse (mean) age female married ptotinc logptotinc aginc wkh fe_cz=__hdfePLACE__ fe_id=__hdfeINDIVID__ logptotinctwfe (count) n_obs_2017_ur=age n_obs_acs_ur=ptotinc n_obs_twfe_ur=__hdfeINDIVID__, by(stfips) wildparse 
		}		
		
		ds stfips *ur, not
		foreach var in `r(varlist)' {				
				  rename `var' mean_`var'
			     }
		gen sample="`occ'; year 2017; age 40 to 55"
		order sample
		save ${temp}/state_`occ'_age40to552017.dta, replace 
		
	}
	
	use ${temp}/state_physicians_allage2017.dta, clear
	append using ${temp}/state_physicians_age40to552017.dta
	append using ${temp}/state_lawyers_allage2017.dta
	append using ${temp}/state_lawyers_age40to552017.dta
	order *ur, last
	replace mean_wkh=. if n_obs_acs_ur<10
	replace mean_income_acs=. if n_obs_acs_ur<10
	replace mean_logptotinctwfe=. if n_obs_twfe_ur<10
	replace mean_fe_id=. if n_obs_twfe_ur<10
	replace mean_fe_cz=. if n_obs_twfe_ur<10
	
	ds *ur sample stfips, not
	drbest_docinc `r(varlist)'
	drbcount n_obs_2017_ur, g(N)
	order sample stfips N
	order *ur, last
	export delimited using "${mypath}/intermediate_csv/desc09-maps_state_level.csv", replace dataf delim(tab)  


********************************************************************************

* CZ level data 
	
	
	foreach occ in physicians lawyers {
			
		*All ages
		
		use if year==2017 using "${clean_data}/panel_`occ'_clean.dta", clear
		format wkh %9.0g	
		
		if "`occ'" == "physicians" {
		gen primary_care=(spec_category_code==1)
		gcollapse (mean) age female married ranked_school top5_school ptotinc logptotinc totw2comp aginc logaginc top1 with25Kbizinc retired_1099ssa wkh income_acs primary_care (count) n_obs_2017_ur=age n_obs_schl_ur=ranked_school n_obs_acs_ur=income_acs, by(cz90) 
		}
		else if "`occ'" == "lawyers" {
		gcollapse (mean) age female married ptotinc logptotinc aginc wkh (count) n_obs_2017_ur=age n_obs_acs_ur=ptotinc, by(cz90) 
		}	
		ds cz90 *ur, not
		foreach var in `r(varlist)' {				
				  rename `var' mean_`var'
			     }
		gen sample="`occ'; year 2017; all ages"
		order sample
		save ${temp}/CZ_`occ'_allage2017.dta, replace 
		
		* Ages 40 to 55 
					
		use if year==2017 & inrange(age,40,55) using "${clean_data}/panel_`occ'_clean.dta", clear
		
		if "`occ'" == "physicians" {
		merge 1:1 id year using "${output}/twfe/geo_twfe/`occ'/IndividFixedEffects_ALL_`occ'.dta", keep(master match) nogen
		}
		else if "`occ'" == "lawyers" {
		merge 1:1 id year using "${output}/twfe/geo_twfe/`occ'/IndividFixedEffects_`occ'.dta", keep(master match) nogen
		}		
		format wkh %9.0g
		gen logptotinctwfe=logptotinc if !missing(__hdfeINDIVID__)
		
		if "`occ'" == "physicians" {
		gen primary_care=(spec_category_code==1)
		gcollapse (mean) age female married ranked_school top5_school ptotinc logptotinc totw2comp aginc logaginc top1 with25Kbizinc retired_1099ssa wkh income_acs primary_care fe_cz=__hdfePLACE__ fe_id=__hdfeINDIVID__ logptotinctwfe (count) n_obs_2017_ur=age n_obs_schl_ur=ranked_school n_obs_acs_ur=income_acs n_obs_twfe_ur=__hdfeINDIVID__, by(cz90) 
		}
		else if "`occ'" == "lawyers" {
		gcollapse (mean) age female married ptotinc logptotinc aginc wkh fe_cz=__hdfePLACE__ fe_id=__hdfeINDIVID__ logptotinctwfe (count) n_obs_2017_ur=age n_obs_acs_ur=ptotinc n_obs_twfe_ur=__hdfeINDIVID__, by(cz90) wildparse 
		}		
		
		ds cz90 *ur, not
		foreach var in `r(varlist)' {				
				  rename `var' mean_`var'
			     }
		gen sample="`occ'; year 2017; age 40 to 55"
		order sample
		save ${temp}/CZ_`occ'_age40to552017.dta, replace 
		
	}
		
	
	use ${temp}/CZ_physicians_allage2017.dta, clear
	append using ${temp}/CZ_physicians_age40to552017.dta
	append using ${temp}/CZ_lawyers_allage2017.dta
	append using ${temp}/CZ_lawyers_age40to552017.dta
	order *ur, last
	  
	 * Measures of income variation at the CZ level
	 
	 foreach var in mean_ptotinc mean_logptotinc  {				
			  gegen sd`var'=sd(`var')
			  gegen cv`var'=cv(`var')
		     }
			
	merge m:1 cz90 using "${intermediate_data}/cz_characteristics/CZ1990Characteristics.dta", keep(master match) nogen keepusing(logpopulation czname)
	gen pop2017=exp(logpopulation)
	rename logpopulation logpop2017
	order cz90 czname pop2017 logpop2017
	
	keep if pop2017>600000
	* As per disclosure rules

	replace mean_wkh=. if n_obs_acs_ur<10
	replace mean_income_acs=. if n_obs_acs_ur<10
	replace mean_logptotinctwfe=. if n_obs_twfe_ur<10
	replace mean_fe_id=. if n_obs_twfe_ur<10
	replace mean_fe_cz=. if n_obs_twfe_ur<10
	
	ds *ur sample cz90 czname, not
	drbest_docinc `r(varlist)'
	drbcount n_obs_2017_ur, g(N)
	order sample cz90 czname N
	order *ur, last
	sort sample cz90
	export delimited using "${mypath}/intermediate_csv/desc09-maps_BigCZ_level.csv", replace dataf delim(tab)  
